# This script runs the regressions required to create Figure S1 of the paper.

library(tidyverse)
library(lfe)
library(ggplot2)
load("dta_single.rdata")
dta_single = dta_single %>%
  filter(year > 2009 & year < 2020)

crime_classes = unique(dta_single$crime_class)
crimes = unique(dta_single$crime)
q = length(crimes)
p = length(crime_classes)
reg_types = c('w_sc_ag', 'w_sc', 'w_int_ag', 'w_int', 'w_shares_ag', 'w_shares', 'w_none', 'w_none_ag')


results_all_int = data.frame(type = NA, weights = NA, agency = NA, estimate = NA, cse = NA, se = NA, N = NA, r2 = NA, adjr2 = NA)

for(i in c(1:q)){
  results_new = data.frame(type = NA, weights = NA, agency = NA, estimate = NA, cse = NA, se = NA, N = NA, r2 = NA, adjr2 = NA)
  ct = crimes[i]
  cat(ct, "\n")
  
  weight_dta = dta_single %>%
    filter(crime == ct) %>%
    group_by(ori9) %>%
    summarize(
      w_sc = mean(subscriberCount, na.rm=T),
      w_int = mean(stats_interactions, na.rm=T),
      w_shares = mean(stats_shareCount, na.rm=T),
      w_none = 1
    ) %>%
    select(ori9, w_sc, w_int, w_shares, w_none)
  
  new_dta = dta_single %>%
    filter(crime == ct) %>%
    group_by(ori9, type) %>%
    mutate(
      b = (black > 0)*1,
      w_ag = 1/n())
  
  new_dta = merge(new_dta, weight_dta, by='ori9', all.x=T)
  
  new_dta = new_dta %>%
    mutate(
      w_sc_ag = w_sc * w_ag,
      w_int_ag = w_int * w_ag,
      w_shares_ag = w_shares * w_ag,
      w_none_ag = w_none * w_ag
    )
  
  for(k in c(1:length(reg_types))){
    reg_type = reg_types[k]
    the_reg = felm(b ~ type | ori9 | 0 | ori9, data=new_dta[!is.na(new_dta[,reg_type]),], w=new_dta[!is.na(new_dta[,reg_type]),reg_type])
    results_new[k,] = c(ct, str_split(reg_type, "_")[[1]][2], grepl('_ag', reg_type)*1, the_reg$beta[1], the_reg$cse[1], the_reg$se[1], the_reg$N, summary(the_reg)$r.squared, summary(the_reg)$adj.r.squared)
  }
  results_all_int = bind_rows(results_all_int, results_new)
}


for(i in c(1:p)){
  results_new = data.frame(type = NA, weights = NA, agency = NA, estimate = NA, cse = NA, se = NA, N = NA, r2 = NA, adjr2 = NA)
  ct = crime_classes[i]
  cat(ct, "\n")
  
  
  weight_dta = dta_single %>%
    filter(crime_class == ct) %>%
    group_by(ori9) %>%
    summarize(
      w_sc = mean(subscriberCount, na.rm=T),
      w_int = mean(stats_interactions, na.rm=T),
      w_shares = mean(stats_shareCount, na.rm=T),
      w_none = 1
    ) %>%
    select(ori9, w_sc, w_int, w_shares, w_none)
  
  new_dta = dta_single %>%
    filter(crime_class == ct) %>%
    group_by(ori9, type) %>%
    mutate(
      b = (black > 0)*1,
      w_ag = 1/n())
  
  new_dta = merge(new_dta, weight_dta, by='ori9', all.x=T)
  
  new_dta = new_dta %>%
    mutate(
      w_sc_ag = w_sc * w_ag,
      w_int_ag = w_int * w_ag,
      w_shares_ag = w_shares * w_ag,
      w_none_ag = w_none * w_ag
    )
  
  for(k in c(1:length(reg_types))){
    reg_type = reg_types[k]
    the_reg = felm(b ~ type | ori9 | 0 | ori9, data=new_dta[!is.na(new_dta[,reg_type]),], w=new_dta[!is.na(new_dta[,reg_type]),reg_type])
    results_new[k,] = c(ct, str_split(reg_type, "_")[[1]][2], grepl('_ag', reg_type)*1, the_reg$beta[1], the_reg$cse[1], the_reg$se[1], the_reg$N, summary(the_reg)$r.squared, summary(the_reg)$adj.r.squared)
  }
  results_all_int = bind_rows(results_all_int, results_new)
}

results_new = data.frame(type = NA, weights = NA, agency = NA, estimate = NA, cse = NA, se = NA, N = NA, r2 = NA, adjr2 = NA)
cat("all", "\n")


weight_dta = dta_single %>%
  group_by(ori9) %>%
  summarize(
    w_sc = mean(subscriberCount, na.rm=T),
    w_int = mean(stats_interactions, na.rm=T),
    w_shares = mean(stats_shareCount, na.rm=T),
    w_none = 1
  ) %>%
  select(ori9, w_sc, w_int, w_shares, w_none)

new_dta = dta_single %>%
  group_by(ori9, type) %>%
  mutate(
    b = (black > 0)*1,
    w_ag = 1/n())

new_dta = merge(new_dta, weight_dta, by='ori9', all.x=T)

new_dta = new_dta %>%
  mutate(
    w_sc_ag = w_sc * w_ag,
    w_int_ag = w_int * w_ag,
    w_shares_ag = w_shares * w_ag,
    w_none_ag = w_none * w_ag
  )

for(k in c(1:length(reg_types))){
  reg_type = reg_types[k]
  the_reg = felm(b ~ type | ori9 | 0 | ori9, data=new_dta[!is.na(new_dta[,reg_type]),], w=new_dta[!is.na(new_dta[,reg_type]),reg_type])
  results_new[k,] = c("all", str_split(reg_type, "_")[[1]][2], grepl('_ag', reg_type)*1, the_reg$beta[1], the_reg$cse[1], the_reg$se[1], the_reg$N, summary(the_reg)$r.squared, summary(the_reg)$adj.r.squared)
}
results_all_int = bind_rows(results_all_int, results_new)

results_all_int = results_all_int[-1,]
results_all_int$estimate = as.numeric(results_all_int$estimate)
results_all_int$cse = as.numeric(results_all_int$cse)
results_all_int$se = as.numeric(results_all_int$se)
results_all_int$N = as.numeric(results_all_int$N)
results_all_int$r2 = as.numeric(results_all_int$r2)
results_all_int$adjr2 = as.numeric(results_all_int$adjr2)
write_csv(results_all_int, 'all_results_alternative_weights.csv')